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Abstract 

Two methods based on Fourier series expansions (a Chan- 
drasekhar functions-based method and a shifted Legendre polyno- 
mials -based method) are used to study analytically the eigenvalue 
problem governing the linear convection problem with an uniform 
internal heat source in a horizontal fluid layer bounded by two 
rigid walls. For each method some theoretical remarks are made. 
Numerical results are given and they are compared with some 
existing ones. Good agrement is found. 
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1 Problem setting 

The effects of the presence in a fluid of an internal heat source have 
been experimentally, numerically and analytically investigated by re- 
searchers in many convection problems [6], [7j, [8],[9]. The investigations 
concerned the effects of the heating and cooling rate. Various conditions 
were imposed on the lower and upper boundaries. The motion in the 
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atmosphere or mantle convection are two among phenomena of natural 
convection induced by internal heat sources. They bifurcate from the 
conduction state as a result of its loss of stability. In spite of their im- 
portance, due to the occurrence of variable coefficients in the nonlinear 
partial differential equations governing the evolution of the perturbations 
around the basic equilibrium, so far these phenomena were treated mostly 
numerically and experimentally. 

Herein a horizontal layer of viscous incompressible fluid with constant 
viscosity and thermal conductivity coefficients v and k is considered |9]. 
In this context, the heat and hydrostatic transfer equations are [9] 
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where rj = const, is the heating rate, 9b, Pb and pb are the potential 
temperature, pressure and density in the basic state, respectively. In the 
fluid, the temperature at all point varies at the same rate as the bound- 
ary temperature, so the problem is characterized by a constant poten- 
tial temperature difference between the lower and the upper boundaries 
A9b = Obo ~ ^Bi- Taking into account ([T]) this leads to the following 
formula for the potential temperature distribution [9J 
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In nondimensional variables the system of equations characterizing 
the problem is 



d\J 

~dt 



-Vp' + AU + Gre% 



{ divU = 0, 



de' 



(4) 



— = (1 - A^z)Uk + Pr-^Ae', 

where U = [u, v, w) is the velocity, 6' and p' are the temperature and 
pressure deviations from the basic state [9J, Gr is the Grashof number. 
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Pr is the Prandtl number and is a nondimensional parameter charac- 
terizing the heating (coohng) rate of the layer. 

The boundaries are assumed rigid and ideal heat conducting, so the 
boundary conditions read 

\J = 9' = 0atz = —- and ^ ~ 2' 

In [0] the numerical investigations concerned the vertical distribution 
of the total heat fluxes and their individual components for small and 
moderate supercritical Rayleigh number in the presence of a uniform heat 
source. 

The eigenvalue problem associated with the equations for a convection 
problem with an uniform internal heat source in a horizontal fluid layer 
bounded by two rigid walls was deduced in [2]. 

Consider the viscous incompressible fluid confined into a periodicity 

rectangular box V : < x < ai, < y < a2, < z < - [3] bounded 

by two rigid horizontal walls. The corresponding eigenvalue problem [2] 
has the form 

(^2 _ a^)2^r _ (PjiaQ = 0, 

{D^ -a^)e + {l- Nz)W = 0. ^ ' 

with the boundary conditions 

W = DW = e = at z = ±-. (7) 

In (EI) the Rayleigh number Ra represents the eigenvalue while (W, 9) 
represents the corresponding eigenvector. The analytical study of this 
stability problem consists in finding the smallest eigenvalue, i.e. the 
critical value of the Rayleigh number at which the convection sets in. 

In [2] the analytical study of the eigenvalue problem ©-© was per- 
formed by means of a method from [1]. First the system ©-(ITj) was 

written in a more convenient independent variable x = z-\ — . Then, two 

methods (one based on Fourier series expansions of the unknown func- 
tions and other a variational one) were used in order to find the smallest 
eigenvalue. Here, the analytical study in also based on Fourier series ex- 
pansions of the unknown functions, but the expansion functions satisfy 
all boundary conditions. 
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Taking into account the form of the boundary conditions two methods 
are used and, for each of them, some analytical remarks on the chosen 
sets of expansion functions are presented. 



2 A method based on Chandrasekhar func- 
tions 

In this method, the unknown function W is expanded upon a com- 
plete set of orthogonal functions that satisfy all boundary conditions 

= DW = at z = =1=2) then, from we find the expres- 
sion of the unknown function B. Replacing these expansions in and 
imposing the condition that the left-hand side of the obtained equation 
to be orthogonal to each function from the expansion set, we obtain an 
algebraic system of equations which leads us to the secular equation, 
yielding the critical value of the Rayleigh number. 

When the normal component of the velocity and its derivative are 

zero at z = —- and ^ ~ 2' classical set of complete orthogonal 
functions that satisfy these conditions are the Chandrasekhar sets of 
functions {Cn}nm, {5'„}„gn[1] 

coshA„,2 cos XnZ 

Cniz) 



coshA„,/2 cosA„/2' 

^ _ sinh(/inz) _ sm{^nz) 
sinh(/i„/2) sin(/i„/2) 

where A„ and the positive roots of the equations tanh ^— j + 

tan ( = and coth — cot = 0. We have 



2/ \2J V2 

0.5 1-0.5 

Cn{z)Cm{z)dz = / Sn{z)Sm{z)dz = 6„ 



0.5 ^-0.5 



By definition, the functions C„ and S'„ and their derivatives vanish 
at X = ±- so the boundary conditions (^^| are satisfied. 
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Let us consider W = WnCn{z). From (I6J2 we obtain the expres- 

n=l 

sion of the unknown function B, 

^ , , , ^ . , , WnCoshA„^(iV^-l) 2A„W„ 
H = /i cosh az + B smh 02; 



(A2 -a2)coshA„/2 (A^ - a2)2 cosh A„/2 

(l-Arz)W„cosA„;2 2A„W„ . 

■ smh AnZ H — r- — — — — sm AnZ, 

(A^ + a^) cos A„/2 (A^ + a^)^ cos A„/2 

wha.eA=„, , and 

(A2-a2)(A2+a2) cosh a/2 

i3 — 



(A2 - a2)2(A2 + a2)2 cosh A„/2 (A2 - a2)(A2 + a^) " 

However, in our case, replacing these expressions in and impos- 
ing the condition that the left-hand side of the obtained equation to be 
orthogonal to Cm, m G N, we obtain an expression in which the physical 
parameter is missing. The mathematical explanation is that the cho- 
sen set of expansion functions introduced an extraparity (inexistent in 
the given problem), leading to the loss of one of the physical parameter, 
in this case the cooling (heating) rate A^. 

Remark. The physical parameter also disappear when the expan- 
sion functions are 5'^, n = 1, 2, .... 

Another explanation could be the fact that we have no physical or 
mathematical reason to assume that W is either even or odd. The general 
form of W will be considered elsewhere. 



3 A method based on shifted Legendre poly- 
nomials 

In order to avoid the loss of A^, we use a different set of orthogonal 
functions, namely a basis of shifted Legendre polynomials (SLP) on [0, 1]. 
Let us modify the system ([H]) by a translation of the variable z, x = 

^ + -, such that the eigenvalue problem becomes 

r {D^ - a^fW - a^RaQ = 0, 

{ {D^ - a2)e + (A^i - Nx)W = 0, 
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with A^i = 1 + — and the boundary conditions 

W = DW = e = at ,T = and 1. (11) 

Starting with the classical Legendre polynomials defined on (—1,1), 
let us introduce the complete sets of expansion functions. We are interes- 
ted in expansion functions that satisfy all boundary conditions. Let 
i^oHO, 1), ^o(0, 1) be two Hilbert spaces [3] 

iJo'(0, 1) = {/I/, /' e L\0, 1), /(O) = /(I) = 0}, 

H'oio, 1) = {/I/, /', r e L'iO, 1), /(O) = /(I) = /'(O) = f (1) = 0} 

and denote by the Legendre polynomials defined on (—1, 1). By means 
of them, we construct the SLP (denoted by us by Qk) on (a, 6), namely 
/ 2x — Qj — b\ 

Qkix) = Lk ( — ; ) . Taking (a, b) = (0, 1) we find that Qk are or- 

V b — a / 

thogonal polynomials on the interval (0, 1), i.e. / QiQjdz = — — -6ij. 

Jo 2^ + 1 

Using the identity [3] 

2i2t + l)Q,iz)=Q[^,iz)-Q[_,iz). (12) 
we define the complete sets of orthogonal functions {0i}i=i,2,... C Hq{0, 1), 

Qi+l ~ Qi-1 



Mz) = / Qiit)dt 



2(2i + 1) 



satisfying boundary conditions <f)i{0) = 0i(l) = at z = and 1 and 

{Ah=i,2,... c//2(o,i), 



= I I Q^+l{t)dtds = ^ 



^0 



Qi+3 " Qi+l Qi+l ~ Qi-1 



(2i + 3)(2i + 5) (2i + l)(2i + 3) 



satisfying boundary conditions /?i(0) = A(l) = Pi{0) = /?j'(l) = at 
z = and 1. ^ ^ 

Remark. We could also work with SLP on (a, b) = ^ ' ~)' 

ever, the choice (a, b) = (0, 1) leads us to simplified numerical evaluations. 
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The system ([U]) can be solved numerically by approximating the so- 
lution (ly, e) by 



iy = J]w,AW, e = J]e,0,(z) (is) 



1=1 i=l 



with Wi and Gj the Fourier coefficients. In this way, the system ([6]) can 
be written in terms of the expansion functions only 

n 

E[W^iD^ - a^f/3i - a^RaQ,(l)i] = 0, 
'Z' (14) 
T.i'diiD' - a')<Pi + (iVi - Nz)Wil3i] = 0. 

i=l 

Multiplying the system by the vector (/3fc, 0^) we obtain the algebraic 
system 

f:[W^/(D2 - a2)2/?,,/?,) - a2i?ae,(0„/?fc)] = 0, 

(15) 

E [6, ({D^ - a2)0„ + iy,iVi(A, 0fc) - iy.iV(^/5,, </>fc)] = 0. 
i=i ^ ^ 

Taking into account the fact that the coefficients W^j, 0j are not all null, 
i.e. the Cramer determinant vanishes, the secular equation has the form 



((D2-a2)2/5„/5fc) -a2i?a(0„/3fc) 

iVi(A, 0fe) - iV(^A, 0fc) ((^2 _ ^2)0,, 0,) 



0. (16) 



The scalar products from ffTUl) are given in the Appendix. 

The system ffTUj) has variable coefficients (functions of x). In this 
case, the following recurrence relation was used for the numerical study 

2^Q* = + + TTT^Qi-l- (17) 

2t + 1 2t + 1 

4 Numerical results 

Taking n = m = 1 we obtained a first approximation of the Rayleigh 
number, which proved to be a good approximation compared to the one 
obtained in [2J. The obtained numerical results are presented in Table 1 



7 



in comparison with the results from ^ . The disadvantage of this method 
is given by the fact that the approximations are hmited by the difficult 
evaluation of the associated matrix for a large number of functions in 
the expansion sets. However, the expressions of the neutral manifolds 
are easy to obtain with this method. For a large number of terms in 
the Fourier series expansions, we must use an algorithm for solving the 
algebraic equation (I16p . For instance in [3], the Arnoldi algorithm is 
used. 



AT 


^2 

a 


Ra — Fourier 


Ra — var.meth. 


Ra — Legendre 





9.711 


1715.079324 


1749.97575 


1749.95727 


1 


9.711 


1711.742588 


1746.804944 


1746.809422 


2 


9.711 


1701.891001 


1737.45025 


1737.450242 


1 


10.0 


1712.257687 


1747.29100 


1747.290998 


4 


10.0 


1664.341789 


1701.62704 


1701.627037 


4 


12.0 


1685.422373 


1723.62407 


1723.624047 


8 


12.0 


1547.460446 


1590.19681 


1590.196769 


9 


12.0 


1508.147637 


1551.72378 


1551.723746 


10 


12.0 


1468.449223 


1512.69203 


1512.691998 


12 


12 


1389.837162 


1434.90396 


1434.903926 


16 


12 


1243.442054 


1288.50149 


1288.501459 


10 


9.0 


1482.527042 


1525.59302 


1525.593072 


11 


9.0 


1446.915467 


1490.55802 


1490.558078 


12 


9.00 


1411.401914 


1455.48233 


1455.482384 



Table 1. Numerical evaluations of the Rayleigh number for various 
values of the parameters and a. 



When the wavenumber is kept constant an increase in the heating 
(cooling) rate parameter leads to a decreasing of the Rayleigh number. 
When = the problem reduces to the particular case of Rayleigh- 
Benard convection and the numerical evaluation lead us to a value similar 
to the classical value for the Rayleigh number, i.e. Ra = 1749.95727 for 
a = 3.117. 
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5 Appendix 

Let us give the expressions of the scalar products occurring in (|T6 
Since in ( fTOl) i the expression ((D^ — a'^)^Pi,Pk) is written as 



let us simplify these products or simply evaluate them. Taking into 
account the definition of the scalar product on L^(0, 1), i.e. (/, (?) = 

fgdz and the boundary conditions satisfied by the expansion func- 



tions, we have 



1 



and 



2i + 3 
1 



i = k, 



2(2i + l)(2i + 3)(2i + 5) 



if i = /c, 



1 



4(2i-l)(2i + l)(2i + 3) 
otherwise 



ifi = k + 2, 



;i8) 



(19) 



Given the fact that = ■ , , 

the product and we get 



(pk+2 



we first evaluated 



2{2i - l){2i + l)(2i + 3) 



■>i, Vk 



4(2i + l)(2i + 3)(2i + 5) 
^ otherwise 



if i = A;, 



if i = )fc - 2, 



(20) 
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Using fpUl) we have 



8{2i - l){2i + l){2i + 3){2i + 5)(2i + 7) 



4(2i + l)(2i + 3)(2i + 5)(2z + 7)(2z + 9) 
1 

16(2i + 3)(2i + 5)(2i + 7)(2i + 9)(2i + 11) 
otherwise 



if z = A;, 



if i = A; - 2, 



if i = /c - 4, 



We also used (l20l) to deduce i.e 

1 



(21) 



A) 



2(2A: + 3) 



</'fc+2) - ('/'i,'Afc)] 



8(2i-l)(2i + l)(2i + 3)(2i + 5) 
3 



if i = A;, 



8{2i - 3){2i - l){2i + l){2i + 3) 
1 

8{2i + l){2i + 3){2i + 5){2i + 7) 



8{2i - 5){2i - 3)(2i - l)(2i + 1) 
otherwise 



if i = /c + 2, 



iii = k-2, 



if i = A; + 4, 



(22) 



Let us remark that {Pi,4>k) = {4>k,Pi)- 

The computations of {D'^cpi, 0^) was simphfied by the expressions of 
the Si functions. We have 



(Qi, Qk 



2i + l 
otherwise 



if z = A, 



(23) 



All the obtained expressions ffT^ - ([21]) are based on the orthogonality 
relationship between the SLP. In deducing the expression below we also 
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used the recurrence relation ffTTl) 



i + 4 



16(2i + 3)(2i + 5)(2i + 7){2i + 9)(2i + 11) 
1 



if z = /c — 5, 



16(2i + 3)(2i + 5)(2i + 7)(2i + 9) 



if i = - 4, 



1 



16(2i + l)(2i + 3)(2i + 5)(2i + 9) 
3 

16(2i + l)(2i + 3)(2i + 5)(2i + 7) 



16(2i - l)(2i + l)(2i + 3)(2i + 5){2i + 7) 
3 



if i = A; — 3, 
if f = A; - 2, 

a i = k — 1, 



16(2i - l)(2i + l)(2i + 3)(2i + 5) 
1 

16(2i - 3)(2i + l)(2i + 3){2i + 5) 
1 



if i = k. 



ifi = k + l, 



i{i = k + 2, 



16(2i - 3)(2i - l){2i + l){2i + 3) 

i + l 

16(2i - 5)(2i - 3){2i - l){2i + l){2i + 3) 

otherwise 



if i = A; + 3 



(24) 



6 Conclusions 



In this paper we performed an analytical study of the eigenvalue prob- 
lem corresponding to a convection problem with uniform internal heat 
source. We pointed out some aspects of the spectral methods that we 
employed concerning the sets of the expansion functions that can be used 
to an analytical study of this problem. As in this case the expansion sets 
of Chandrasekhar functions introduced an extraparity they were not ap- 
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propriate. However, for some other problems ^ their use proved to be 
successful. The method based on SLP lead to good numerical approxi- 
mations. All numerical results obtained with this method are compared 
with the existing ones. The effect of the heating (cooling) rate on the 
values of the Rayleigh number is pointed out. 
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